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Abstract 

A finite element program is presented to simulate the process of packing and coiling elastic wires in 
two- and three-dimensional confining cavities. The wire is represented by third order beam elements 
and embedded into a corotational formulation to capture the geometric nonlinearity resulting from large 
rotations and deformations. The hyperbolic equations of motion are integrated in time using two different 
integration methods from the Newmark family: an implicit iterative Newton-Raphson line search solver, 
and an explicit predictor-corrector scheme, both with adaptive time stepping. These two approaches 
reveal fundamentally different suitability for the problem of strongly self-interacting bodies found in 
densely packed cavities. Generalizing the spherical confinement symmetry investigated in recent studies, 
the packing of a wire in hard ellipsoidal cavities is simulated in the frictionless elastic limit. Evidence is 
given that packings in oblate spheroids and scalene ellipsoids are energetically preferred to spheres. 

Keywords: Beam, Bending, Corotational formulation, Finite element. Large deformation, Nonlinear, 
Wire, Contact 



1. Introduction 

The dense packing and crumpling behavior of long, slender objects subject to spatial confinement is 
a recurring topic in various fields in natural sciences. Biophysicists that find long DNA chains densely 
packed in viral capsids (e.g., Kindt ct al., 2001; Grayson and Molineux, 2007; Roos et al., 2007; Petrov 
and Harvey, 2008) or surgeons that treat saccular aneurysms in brain arteries by means of endovascular 
coiling (Guglielmi et al., 1991a, b) are just two out of many examples why the coiling of thin rods inside of 
a confined space has attracted increasing interest in recent research. In particular the minimally invasive 
surgical treatment of saccular aneurysms, where a platinum wire is pushed through a small opening of 
the sphere-like bulge, calls for a deeper understanding of the emerging morphologies, since high packing 
densities favor a better long-term stability of the embolization (Tamatani ct al., 2002). 

In the past decade, the crumpling of thin wires has been studied in two dimensions experimentally by 
Donato ct al. (2002, 2003, 2006); Donato and Gomes (2007); Gomes et al. (2008b) followed by numerical 
simulations by Stoop et al. (2008), from a primarily physical perspective. It wasn't until very recently that 
experiments were extended to three dimensions (Gomes et al., 2008a; Stoop et al., 2011). A numerical 
study of the three-dimensional case by Stoop ct al. (2011) has unveiled a range of interesting morphological 
phases largely dependent on boundary conditions and internal twist rather than the stiffness or amount 
of friction of the wire. Stoop et al. used a discrete element model to represent the wire dynamics. 

In this paper, a higher-order representation for elastic wires packed inside of two- or three-dimensional 
cavities is presented, using a finite element scheme. The finite element model is more flexible toward adap- 
tive mesh refinement and adaptive damping, and further amenities include a less involved incorporation 
of plasticity effects, as well as guaranteed convergence. Using the presented model, we generalize previous 
studies on spherical cavity shapes to ellipsoidal ones. The results are indicating that wires coiled inside of 
flexible cavities (e.g. biomaterial antra) will tend to assume non-spherical bulk shapes, differing to some 
degree from the emerging morphologies recently found in perfectly spherical confinement. 

The paper is organized as follows: Section 2 presents the construction of a third order beam finite 
element model representing an isotropic wire that is pushed into (or pulled out of) a hard ellipsoidal 
cavity. Both theoretical and technical aspects of the geometrical nonlinearity arising from large rotations 
and deformations are addressed and amended by the incorporation of wire-cavity contact and wire-wire 
contact. The main components of the strain energy are expressed in terms of the finite wire elements. 
Section 3 addresses the integration of the hyperbolic equations of motion in time, using an implicit and 
an explicit scheme of the Newmark family (Newmark, 1959). The presented beam theory is verified by 
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comparison to literature in Section 4, where we also argue that a purely implicit treatment of wire-wire 
contacts within the finite clement framework is unfeasible in the context of dense packings. Section 5 
briefly presents new simulation results for a selection of different ellipsoidal isochoric and isoareal cavity 
shapes, unveiling their respective energy levels. In Section 6, the results are summarized and an outlook 
to future extensions of the present work is given. 



2. Finite element model 



2.1. Equations of motion 

The wire is represented by a one-dimensional regular mesh with N elements {e}, i.e. N + 1 nodes {n}, 
with uniform spacing h. The undeformed mesh is positioned at x e aligned to the global x-axis 

in the simulations. Each mesh node carries six degrees of freedom (DOFs): three Cartesian displacements 
tt„, Vn, Wn in y-, and z-direction, respectively, and three Euler angles tpn, dn, ipn for rotations around 
these respective axes, giving each node an orientation. The center line of the wire is interpolated along 
each element using cubic Hermite shape functions. The motion of the wire is described by the nonlinear 
system of hyperbolic differential equations 



Mii+Cu + fi„t(u) =fe,t(u). 



(1) 



Here, M is the wire's mass matrix, C is a damping matrix and fint(u) represents the solution-dependent 
internal elastic forces. Each of the matrices is an element of K6(w-hi)x6(jv-hi)^ -pj^g gig^g^j solution vector 
u = u{t) S k6(A'+i) consists of the six DOFs for each node at time t. All external forces and torques 
acting on the wire at time t enter the solution-dependent right-hand side vector fext(u). This includes 
self-contact forces between wire elements and forces due to repulsion from the cavity. The conventional 
DOF ordering u = [i 



^1 , ■ 



u^] is adopted in the following, with nodal contributions 



0,1, ...,7V. 



(2) 



To calculate local element contributions, as is common in the finite element method, we employ the local 
element ordering 

Ue = [ui,Vi,Wi,ipi,9i,lpi,U2,V2,W2,(p2,02,'lp2]'^ G M^^, (3) 

where the subscript 1 refers to the element's left node and 2 its right counterpart. In order to avoid the 
singularities occurring in Eulerian representations of large rotations, auxiliary nodal unit quaternions q„ 
are used in addition. They are initialized according to 
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(5) 



storing the orientation of each mesh node n in space. Given a small increment Au„ to (2), be it obtained 
in whatever way, the quaternion can be updated by applying the non-commutative Hamilton product 



q« = Ago 90 - (Aq)'^q + i {qo Aq + Ago q -I- Aq x q). 



(6) 



where Aq^j is built using Eq. (4) on the angular increments A(p„, A0„, A-f/jn from Au„. The purpose in 
defining triads as in Eq. (5) is elucidated in Section 2.3, where the nodes may be arbitrarily oriented to 
allow for very large deformations. 

For simplicity, the mass matrix M is diagonalized by means of direct mass lumping (Tong et al., 
1971): 

M„ = diag(m„,TO„,m„, 

Jn: Jji: Jn) (7) 



with nodal masses m„ = Ahp^ where A 



is the cross-section area of the tubular wire and p denotes its 



mass density. For the nodal moments of inertia J„, the point masses m„ are uniformly distributed within 
balls of the same radius as the wire, resulting in J„ = |m„r^. For damping, one may use the Rayleigh 



ansatz C = aM -t- /3K, a,/3 > 0. In the present simulations, this is simplified to C = cl, 



where 



c > is a scalar viscous damping coefficient acting equally on all translational and angular velocities, 
and Ife denotes the kxk identity matrix. 
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2.2. Three-dimensional beam theory 

The bending deformation of the tubular wire is modeled by Reddy's simplified third-order beam theory 
(RBT) (Reddy et al., 1997; Reddy, 1997) that includes a quadratic transverse shear stress distribution, 
contrary to previously published discrete element simulations (Stoop et al., 2008, 2011). RBT completely 
alleviates the problem of shear locking (Tessler and Dong, 1981; Prathap and Bhashyam, 1982) as it is 
known in Timoshenko beam theory (Timoshenko, 1921, 1922), without complicating the stiffness matrix. 
Assuming independent bending in both the xy and xz planes, in particular vanishing mixed-bending 



Cauchy stresses cr„ 



0, the two-dimensional linear RBT element stiffness matrix for constant 



beam material and geometry can be extended to 3D as follows: 
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is the parameter determining the order of the theory: 
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In RBT, the stiffness coefficients 
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(i) . 

are needed, in which r^^ is the half thickness of the beam in z direction. lyy is the i-th moment of inertia 
of the wire about the y-axis, and is calculated as 



(13) 



(14) 



Hi) ^ 

yy 



z' AA. 



(15) 



Relations (11-15) hold analogously with y and z interchanged. E denotes the Young's modulus of the 
isotropic wire, G — E/2{1 + ly) the elastic shear modulus, and the Poisson ratio. The longitudinal and 
torsional spring constants are given by 

EA , GJ 
X' 



ku — 



k(p — 



(16) 



where the polar moment of inertia reads J 
with circular cross-section and radius r 



lyy + Izz by the perpendicular axis theorem. For a wire 
- /(O^ which imphes that 



ry = r^, one can write lyy 



kv — k^j ^ kip — ke^ ky^p — 
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(17) 
(18) 



such that fi^jy = 0,^^ = x§g(l + i^) (x) in RBT. Note that in Eq. (8) the rotational DOFs in the two 



bending directions have opposite sign convention such that both "slopes" 6 and (j) are equally oriented 
as their respective conventional Euler angles. The longitudinal displacement u and the twist angle Lp are 
modeled with Hooke's law. 
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2.3. Corotational formulation 

In order to comprise large rotations and deflections in a finite element simulation, geometric nonlin- 
earity must be introduced to the theory, which can be expressed by the fact that the dependence of the 
internal forces fint on u is not a mere linear product with a constant stiffness matrix K. To this end, three 
prominent approaches have been developed, refined, mixed, and broadly applied over decades. The total 
Lagrangian formulation (TL) expresses the sought solution with respect to an initial configuration. The 
updated Lagrangian formulation (UL) can be shown to be effectively equivalent to TL (Shabana, 2008), 
but makes use of a different representation, where the nodal displacements and rotations are described 
incrementally. The corotational formulation (CR) introduces a local reference frame for each element 
that continuously moves and rotates with the element, within which the theory can be considered linear. 
CR has been acclaimed superiority over TL and UL in accuracy (Mattiasson ct al., 1986; Urthaler and 
Roddy, 2005), performance (Belytschko and Hsich, 1973; Rankin and Brogan, 1986), simplicity (Teh and 
Clarke, 1998) and generality (Rankin and Brogan, 1986; Crisfield, 1990), and it is utilized in the following. 

The corotational formulation is rooted in in the works of Wempner (1969), Belytschko and Hsich 
(1973), and Oran (1973a, b). Assuming large displacements and rotations but small strains, the idea is to 
corotate a dedicated coordinate frame with each beam element, exactly passing through the nodes. The 
12 global DOFs per element are reduced by its five rigid body DOFs to seven local deformation DOFs 
when expressed with respect to this corotated frame. Crisfield (1990, 1997) was the first to provide a 
beam-theory-independent, three-dimensional CR formulation with consistent treatment of the tangent 
stiffness matrix required by the Newton-Raphson method. Details of the algorithm can be found there, 
here, only the application to RBT wire elements is briefly discussed. 

Given the global nodal positions p„ = x„ + [u„, u„, and orientations T„, each beam element e 
is assigned a unit triad Eg = [61,62,63] in the CR setup, called the corotated frame, in such a way that 
it represents the orientation of e. In particular, ei connects the two nodes: 



ei 



P2 - Pi 
|P2 - Pll 



(19) 



The global element DOFs Ug are then expressed with respect to Eg . In the following, all variables carrying 
a hat refer to this local element frame. By construction, we have vi — V2 — wi ^ W2 ^ 0. Additionally, the 
longitudinal displacement variables ui, U2 can be reduced to a single longitudinal spring-like deformation 
'^12 = |P2 — Pi I — ^- It turns out that if 62, 63 are constructed from the nodal quaternions according to 
Crisfield, the remaining six local rotational DOFs with respect to Eg can be computed as 



(^„=sin 1 (i(tJe3-tTe2)) -(^°, 
^„=sin-i(i(eTti-tTei))-^Jl, 
7^„-sin-i(i(eTti-tTei))-V>o, 



(20) 



for both nodes n adjacent to element e, each using their respective triad T„ = [ti,t2,t3]. ip'^ , 0'^ , ip'^ 
allow to supply the wire with intrinsic curvature or twist, i.e. curvature or twist occurring without 
external loads or moments. For instance, a uniform intrinsic curvature k in the xz plane is achieved 
by setting —9i — 62 — hK/2 in each element frame. This defines the seven local deformation DOFs 
Ug = [i^i, ^1, i/'i, U12, <^2, ^2, from which the rigid body motion has been completely detached, and 
which are shown together with the triads in Fig. 1. The corresponding 7x7 local element stiffness matrix 
acting on them is easily devised from Eq. (8) using the definition of the corotated frame: 
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(21) 



While the local element stiffness matrix remains constant throughout the simulation, nonlinearity enters 
the equations of motion when Kg is transformed into the desired 12x12 "nonlinear" global element 
stiffness matrix Kg(ug) via a solution-dependent global-to-local transformation matrix rg(ug) € 



p7xl2. 



Kg(ug) =Fj(ug)Kgrg(ug). 



(22) 



Equivalently, one can write 



fint,e = Fj(Ue)KeUg 



(23) 
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for the internal force vector. Note that Ug contains the intrinsic angles ifP,9^,ip'^ via Eq. (21), such that 
fint,e may be non-zero even if the element solution in global coordinates, Ug, vanishes. The construction 
of Fe is done row-wise and comprises somewhat lengthy matrix calculus involving the three triads. 




Figure 1: Visualization of tlie corotated wire element: In the initial straight configuration (a), each element carries canonical 
triads Ee = Ti = T2 = I3. After deformation (b), the seven local element DOFs (bearing a hat) are found with the help of 
the triads Ti (red) and T2 (blue) storing the nodal orientations, and the corotated element frame Ee (magenta) encoding 
the element orientation. 



The most significant part of Crisfield's work is the derivation of a tangent stiffness matrix 

Kt,e(Ug) = -^Ke(Ug)Ue = Ke(Ue) + K,,e (24) 

that is consistent with the rest of the formalism, which is needed for implicit integration of the equations 
of motion in time with the Newton-Raphson method. The geometric stiffness matrix K^-.e is an intricate 
sum of submatrices expressed in terms of the triads, the internal forces (23), the local DOFs Ug and the 
corotational transformation matrix Fg. Details are readily available in Crisfield (1990, 1997). 

This setup allows to numerically integrate the solution vector u in time. In each time step, and 
in each iteration of the nonlinear solver routine, for each element, the current corotated frame Eg, the 
seven local DOFs Ug, and the transformation matrix Fe are computed. The internal force (23) is then 
calculated, and if required by the integration method, the tangent stiffness matrix is built using Eqs. (22) 
and (24). After a new solution vector has been found with the applied integration scheme, the nodal 
quaternions are updated by applying the quaternion product (6), where the Aq„ are built using Eq. (4) 
on the vectorial difference between new and old solution. The quaternions need to be renormalized on a 
regular basis. 

2.4^. Strain energy 

For Euler-BernouUi theory, it can easily be shown by integrating the cubic Hermite splines over the 
corotated elements, that the potential energy due to bending can be expressed in terms of the local DOFs 
by 

c/b = ^ ^ {ol + eA + el + + ^^^2 + ^i) , (25) 



in which subscript 1 and 2 again refer to the left and right node of element e. In RBT, however, the 
Hermitian shape functions are modified to account for the shear strain approximation, and Eq. (25) is 
valid for relatively thin wires only. In the packing simulations presented in Section 5, bending shear 
effects turn out to be small enough to use Eq. (25) as a very close approximation. Independent of the 
employed bending theory, the axial elastic strain energy of the wire is straightforwardly obtained from 
Hooke's law in the corotated frames: 

t^^ = |^E"?2. (26) 

e 

Analogously, the torsional potential energy reads 

U^-^Y.^V2-^if. (27) 
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2.5. Wire-cavity contact 

Spatial confinement is introduced by a rigid ellipsoidal container centered at the origin, with a single 
small opening through which the wire is pushed inside. Hertzian contact mechanics is used to model 
the repulsive forces exerted on the wire nodes by the cavity. In the event that the sphere with radius r 
around a node n indents the cavity wall by a depth £> > 0, the normal force f± — E*hD-K/A (Johnson, 
1985) is applied to the three translational DOFs of the penetrating node, i.e. the nodal contribution to 
fext is 

fcav,ri = — /_L n„ G M? . (28) 

Here, n„ is the outer normal vector of the cavity surface at the contact point, and E* mixes the material 
parameters (£', v) of the rod with those of the cavity, (-Bcav; J^cav): according to 

1 \ — v"^ 1 — v]:„^, 
— = \ (29) 

E* E i?cav 

The normal force is that of a Hertzian contact between two parallel cylinders, which serves as a reasonable 
rough estimate for the real force exerted by a large, hard cavity to a cylinder. Alternative models such 
as the Hertzian force between a sphere and a half-space, where the force is not linearly proportional to 
D, may be used instead if desired. The exact kind of the applied force law turns out to be of very minor 
importance for the present packing model. 

The calculation of the shortest distance between an exact ellipsoid and a point in space is known to 
be equivalent to the problem of finding the roots of a sixth order polynomial (Hart, 1994), for which no 
analytical solution is known. Thus, the indentation depth D can only be found numerically, and so can the 
Jacobian, which needs to be specified in the Newton-Raphson scheme for truly implicit integration in time. 
Particularly concerning the Jacobian, this is very inconvenient. Instead, a closed-form approximation is 
developed here. Consider an ellipsoidal cavity given by 

For reasonably shaped ellipsoids not too far away from a sphere, i.e. for ~ Ry ~ Rz ^ r, a cavity 
contact is assumed if A > in the implicit formula of an effective ellipsoid 

given the nodal position p„ = [px,Py,Pzf^ ■ Since the Hertzian force outgrows the internal forces of the 
wire even for small indentations D when compared to the length scales of the wire, its magnitude is far 
less important than its direction. Therefore, D can legitimately be approximated using 

(R-r + Df^pl+pl+pl^(R-rf{l + A), (32) 

where R — {R^ + Ry + R^) /3 averages the ellipsoidal radii. Taking the square root yields the sought 
approximation 

D ^ (R-r){VTTA^l). (33) 

A consistent approximation of the outer surface normal vector assuming that A ^ 1 is found by the 
normalized gradient of the effective ellipsoid: 

Vgeff(Pn) _. Pn^ . . 

iv£;eff(p„)r- ip„r 

The effective direction p„ can be compactly written as p„ = 2 R^^p„ with a diagonal matrix 

R = diag(i?:r -r,Ry- r, R^ - r) £ M?'^'^ . (35) 

The Jacobian Jcav.n G M?'^'^ of the nodal force vector is easily assembled using the standard rules of 
differential calculus: 

, i9fcav,n , / 9n„ n„ ai) \ 

Jcav,n = — -J_L ^ H — ^ , (,ODj 

dp„ \OPn D dpnj 

^ « — f I - — 
dpn ^ \Pn\ \ ^ \Pn 

dD R-r ,^ 

■TT— ~ , p^. 38 

dpn 2\/rTA " 



PnPl n ^ (37) 
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The impact of the approximations in Eqs. (33) and (34) is twofold. First, the magnitude of the 
Hertzian contact force is subject to an error that grows with the disparity of R^, Ry, Rz- Since the 
Hertzian contact model used here is itself only a crude approximation, this has no negative consequence 
in practice. The only purpose is to keep the wire inside the confining container. Second, the effective 
ellipsoid has slightly deformed aspect ratios compared to the aimed cavity. This defect is also minor for 
Rx ~ Ry ~ Rz and vanishes for arbitrary R^, Ry, Rz in the slender wire limit r — > 0. All approximations 
above become exact under spherical symmetry (R^ ~ Ry — Rz)- 

The axial symmetry of the straight wire is broken by imposing a small initial random deflection 
perpendicular to the wire on two nodes. The mesh is dynamically extended by additional elements as the 
wire enters the cavity. Boundary conditions are imposed by constraining nodes outside the cavity to the 
X-axis, penalizing rotations in the outermost node (no dissipation of torsional energy out of the cavity), 
and by applying a constant insertion velocity un = vin- The packing density is computed using 

'^=^nr i'^^D, 0^ . „ p „ in 3D, L = J2ii^i2 + h). (39) 
■kRxRz 4:/3nRxRyRz ^ 

2.6. Wire-wire contact 

Self-interaction of the wire happens among pairs of wire elements (ei, 62), and the resulting forces act 
on the translational DOFs of the four involved nodes ni,n2 of element ei and 713,714 of element 62. For 
simplicity, all wire elements are treated as if they were sphero-cylindrical in shape. Since the number 
Nc of element pairs in contact grows like the square of the packing density in the uncorrelated thin rod 
limit in 3D (Philipse, 1996) (see also inset of Fig. 7b), an efficient treatment of these pairs is crucial for 
performance. Calculating the shortest distance between any two elements would lead to 0{N^) distance 
calculations, but only a small subset of them will actually be in contact. We therefore separate the 
problem into two steps. 

First, possible contact pairs are found efficiently by using linked cell lists (Quentrcc and Brot, 1973). 
The linked cell algorithm assigns every element to a cubic cell, such that only elements assigned to the 
same or adjacent cells have to be considered as possible contact candidates, reducing the complexity to 
0{N). The minimum cell size guaranteeing that any contact between two wire elements happens within 
neighboring cells is I = 2{h + r). To avoid having to rebuild the linked cells after each iteration step, a 
small additional margin can be added to the cell size, and the cells are rebuilt only after the wire has 
moved farther than half of that margin. 

In the second step, the shortest distance between every contact candidate is calculated using the 
quadric algorithm by Sunday (2001), which provides an efficient implementation based on a note by 
Ebcrly (1999, 2000), with the modification that the center line segments are normalized to unit length to 
fix the faulty detection of short parallel segments in Sunday's algorithm. The algorithm finds si, S2 £ [0, 1] 
such that the closest points of approach on the wire center lines read 

Cl = P„i -I- Si(p„2 - p„J, C2 = p„3 S2(Pn4 - Pns)- (40) 

The mutual indentation of the elements is then given by its depth D = 2r — |Ac| and direction n = 
Ac/|Ac|, where Ac = Ci — C2. 

Similar to cavity contacts, the spatial self-interaction of the wire is modeled by exchange of Hertzian 
contact forces that enter the equations of motion (1) through the right-hand side foxt- No moments are 
transferred, i.e. all contacts are considered frictionless. The same normal force f± as for cavity contacts 
is used, which is motivated by the observation that most self-contacts occur between almost parallel wire 
elements. In favor of a compact and general notation, we condense the four involved nodal positions into 
a single vector 

P = [pI,,P^,,P^3'PnJ^e]R'', (41) 

and apply the same ordering for the resulting force vector fscif,ei,e2 G ^^'^ and its Jacobian Jsoif,ei,e2 G 
1^12x12 -j-j^g following. The Hertzian normal force is distributed to the four nodes by setting 

fsolf,ei,e2 =/-l(w«)n) (42) 

with linearly interpolated nodal weights w = i[l — si,si,— (1 — S2),— 52]""" (see Fig. 2). (E) denotes the 
Kronecker product. In other words, the closer a node lies to the contact point between the two elements, 
the larger the force it experiences, pushing it away from the other element. All possible types of element- 
element contact are covered by this unified description, be it between coplanar elements, or at any angle. 
It also inherently handles the cases where the contact point lies on a spherical cap of at least one of the 
two elements, instead of on their lateral surfaces. Three of the numerous possible contact scenarios are 
visualized in Fig. 2. 
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Figure 2: (a): Wire self-contact, with closest points of approach (blue). The repulsion force (red) is distributed to the four 
involved nodes. (b)-|-(c): A small shift in one of the elements can cause the contact parameters si, S2 to jump, leading to 
singularity in the force Jacobian. 



For a Newton-Raphson solver fully accounting for changes in contact forces during each iteration, the 
derivative of Eq. (42) with respect to a change of the nodal positions needs to be specified. By the rules 
of differential calculus, it reads 

J..„,„,« = *Sil^ .(wanj^t/xf^Sn + wa^^e R'"", (43) 

^ = (wT»l3) + (p,„ - p,„)^ _ (P". -P".)|| ^ («) 

From Eqs. (43-46), the Jacobian matrix is fully defined, given that dsi/dp and 832/ dp can be calculated. 
However, these are singular if the two wire elements are parallel, since si and S2 may jump from one node 
to the other (see Fig. 2), causing convergence difficulties in the Newton-Raphson method. The exchange 
of forces between contacting elements is therefore further simplified in the case of implicit integration 
in time by distributing the normal force equally to all four nodes, i.e. w = [ji |:j ^ j] ? a-nd by 
discarding dsi/dp = 882/ dp = 0. 



3. Integration in time 

3.1. Newmark methods 

Newmark's family of integration methods (Newmark, f959) is widely used for solving hyperbolic 
equilibrium equations in structural dynamics (Zeng et al., 1992; Rcddy, 2004). ft is particularly well suited 
in our case because very simple error estimators are available for the predictor-corrector counterparts of 
the implicit schemes. Although its high precision may not strictly be necessary to account for the dynamic 
relaxation induced by our choice of lumped masses and damping, it is assisting in efficiently and stably 
resolving the steep contact forces. For fixed scheme parameters (3 and 7, the solution vector and its time 
derivatives are discretized and integrated according to 

u{t + At) w ut+At =ut + Atut + ((1 - 2/3)ut + 2/3 iit+At) , (47) 

u{t + At) w Ut+At ^ut + At ((f - 7)ut + 7 iit+At) • (48) 
At denotes the finite time step, and the nodal accelerations at time t are given by 

ii(t) « lit = M-i(fe,t(ut) - fi„t(ut) - Ciit). (49) 

3.2. Implicit integration 

In the case that > 0, the system of equations of motion (I) reduces to the system of nonlinear 
algebraic equations (Reddy, 2004) 

rt,t+At = Kt+AtUt+At - it,t+At = (50) 
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with 



Kt+At - aoM + aiC, (51) 

ft,t+At = foxt(Ut+At) - fint(Ut+At) + Ma^ + Cbj, (52) 

at = ao Ut + a2 lit + 03 lit, (53) 

bt = aiUt +a4Ut + a5Ut, (54) 



and Newmark coefRcients 



""'^jikif' "'^jp-^' = '^^ = ^K$"'^- ^^^^ 



In the present huplementation, the unconditionally stable constant-average acceleration method, that is 
obtained by setting /3 = 1/4 and 7 = 1/2, is chosen. Our Newton-Raphson (NR) line search algorithm 
solves Eq. (50) by taking full account of updated contact forces during the iterations. This requires 
calculating the residual vector and the Jacobian matrix based on the current iterated guess for Ut+At- 
The residual is given by Eq. (50), where fext is the assembly of nodal contributions from cavity contacts 
(28) and self-contacts (42), while fint is the assembly of corotated elastic element forces (23). The Jacobian 
is built according to 

Jt,t+At = Kt + flo M + fli C - Jcxt, Jcxt = Jcav + Jsolf (56) 

in which Kt is the assembly of element tangent stiffness matrices (24), Jcav the assembly of cavity 
contact Jacobians (36), and Jsoif the assembly of self-contact Jacobians (43). Upon convergence to the 
new solution Uf+At, the corresponding velocity and acceleration vectors are calculated using 

ii-t+At = ao(ut+At - ut) - a2 lit - 03 lit, (57) 
lit-HAt = lit + At ((1 - 7)ut + 7 iit+At) ■ (58) 

In a simple attempt to optimize performance and to avoid divergence dead ends in energetically difficult 
configurations, we half the time increment At as long as more than a predefined number of line search 
iterations are required for convergence, or double it provided that the number of iterations remains small. 

Note that an alternative approach, where contact forces are kept constant during the NR iterations, 
is conceivable. This would negate the need to compute Jext, which would just drop out of the formalism, 
greatly accelerating convergence within a single iteration. See Section 4.2 for the reason why this is 
unfeasible in the present context. 

3.3. Explicit integration 

For explicit integration in time, the constant-average acceleration method is implemented in combina- 
tion with an adaptive time stepping predictor-corrector (APC) algorithm based on an a posteriori local 
error estimator by Zienkiewicz and Xie (1991); Zeng et al. (1992). The prediction step is obtained from 
the explicit part of Eqs. (47) and (48): 

(At)"^ 

<+A* = ut + Atiit + ^^(l-2/3)ut, (59) 
uf+At = ut + At(l-7)ut, (60) 



u 



= M-l(fe,t(uP+^t) - fint(uP+^t) - CliP^^,). (61) 

Owing to lumped masses, M^^(-) is a mere component-wise vector multiplication. Based on the predicted 
acceleration, the displacements and velocities are corrected by the previously omitted part: 

^t+At^^^t+At + m'P^t+Av (62) 
u?+ At = At + 7 uf+ At • (63) 
A simple estimator for the relative local error made by performing such a step is then given by 



Vt+At = 



(AtV 1 

^|!uP_,^,-Ut|loo, /3^^, (64) 



where || • ||oo denotes the maximum norm, which is chosen here as a conservative metric. Urcf is a 
characteristic reference length for normalization, such as the cavity size itrcf = R for instance. Ideally, 
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the tradeoff between large time step and small error is dealt with in such a way that the local error 
is approximately constant over time. Zienkiewicz and Xie found that this can be efficiently achieved 
by applying the following adaptive time stepping rules, given a desired target value fj for the relative 
local error, a maximum tolerance jymax > J?, and a lower bound rymin < V above which the time step is 
considered large enough: 

• If T]t+At G [?/min, '7inax], acccpt the time step without modifying At. 

1/3 

• If rjt+At < Vmin, accept the time step and continue with At := (^/jyt+At) At. 

1/3 

• If rjt+At > VmaK, reject the time step and repeat with At := {fj/r]t+At) At. 

4. Verification and performance 

The implementation of the finite element program has been realized with the aid of the libMesh 
(Kirk et al., 2006) and PETSc (Balay et al., 1997) C++ libraries. The applied nonlinear solver is a cubic 
Newton-Raphson backtracking line search algorithm by Dennis and Schnabel (1983). 

4-.1. Test example 

Verifying and benchmarking nonlinear beam implementations can be done by comparison to various 
published static test examples (e.g.. Bathe and Bolourchi, 1979). One of the most prominent ones is 
reiterated here: The 45 degree bend scenario is recommended by the National Agency for Finite Element 
Methods and Standards (NAFEMS, 3DNLG-5) (Abaqus). It involves a cantilever beam with square cross- 
section and non-zero intrinsic rotation variables \0^\ = h/2R. The thereby obtained eighth of a circle with 
radius R = 100 is subjected to transverse end loads Q leading to a three-dimensional nonlinear response. 
We have solved it implicitly and fully static, i.e. with p = c = 0, and depicted it in Fig. 3. The resulting 
tip positions are listed in Tab. 1 for E = 10^, i/ = 0. The reported difference between corotated EBT 
(CR-EBT) and Crisfield (1990) appears to be due to the significantly tightened convergence criterion 
||r||2 £ 10^^ adopted here. Since the order of the theory has almost no effect on such a thin beam, results 
for thickness 10 are also given, revealing the improved precision of corotated RBT over EBT at large 
shear stresses. 



Figure 3: 45 degree bend consisting of 8 elements with square cross-section A, subject to transverse end loads Q. 



4-. 2. Implicit vs. explicit integration 

A general advantage of implicit integration in time over explicit schemes is the ability to cope with 
large time steps without losing numerical stability. In many problems including the present strong spatial 
confinement, however, this benefit is compromised by significantly more time-consuming computations for 
a single time step. The denser a cavity is packed with a wire, the more complex and rough the potential 
energy landscape gets, with fatal consequences for the convergence behavior of an implicit integration 
method. As the NR line search solver needs to find a path through energy valleys getting narrower and 
narrower, there is an increasing chance of starting far enough from the sought energy minimum that 
no path is found any more. In such situations, adaptive time step control is indispensable for avoiding 
divergence dead ends and to ensure the found local minimum corresponds to the physical solution. Fig. 4 
visualizes the implications in practice. A two-dimensional circular cavity with effective size R/r — 50 is 
packed with a frictionless elastic wire, at two different insertion velocities. The "spiral" morphology (Stoop 
et al., 2008) evolving from a slow insertion can be simulated to large packing densities even implicitly at 
relatively large constant time steps At (Fig. 4a). Using the same step size but an insertion speed high 
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CR-RBT (present) 
CR-EBT (present) 


58.25 
58.38 


41.49 
41.22 


22.03 
22.09 


51.54 
51.70 


49.98 
49.67 


18.26 
18.32 


46.35 
46.54 


55.09 
54.75 


15.43 
15.48 



Table 1: Comparison of tip positions of a clamped 45 degree bend at different end loads and beam thicknesses. Initially 
(0 = 0), it is (70.71, 0.0, 29.29). 



enough to introduce dynamic effects from the insertion, the rough valley structure emerging in the local 
energy landscape inhibits the line search solver from converging as soon as wire-wire contacts reach a 
certain degree of complexity (Fig. 4b) . Convergence here means meeting one of the convergence criteria 
within 1000 nonlinear iterations and at line search step sizes A > Amin = 10^^^. The used convergence 
criteria include ||r||2 < 10~^°, ||r||2 < 10^^°||ro||2, and ||Au||2 < , where Tq is the initial residual 

vector, and Au is the proposed NR increment. Simulating the "classic" morphology (Stoop et al., 2008) 
to large densities in reasonable CPU time is within the realm of possibility for the fully updated implicit 
solver only with adaptive time stepping (Fig. 4c). 




Figure 4: (a, 0.8): Low insertion speed, "spiral" morphology. Implicit integration converges even at large constant At. 
(b, 4> Si 0.3): High insertion speed, constant At. Implicit integration convergence failure sets in at very few self-contacts, 
(c, 0.6): High insertion speed, adaptive time stepping. Convergence is retained throughout the simulation. The color 
encodes the wire- wire contact energy from zero (blue) to high (red). For better visibility, the three energy scales are distinct. 

Given the steep energy landscape inherently arising from strongly confined wire packings, it may seem 
tempting to artificially smoothen it by treating contact forces constant during a single NR solver sweep. 
That is, not updating foxt during solving, and hence setting Jcxt = 0, as indicated at the end of Section 
3.2. While this simplifies the formalism and greatly accelerates convergence within each solver sweep, the 
number of such iterations, that is needed per time step until the residual meets one of the convergence 
criteria, is large enough to effectively abolish the gain from faster convergence. Reducing this number 
of iterations to an expedient range requires abridging the time step to unbearably small values. For the 
sum of these observations, we have desisted from such measures, hence keeping the contact forces fully 
updated during all solver iterations. 

Another inherent, structural reason why an implicit integration in time is in general not feasible 
is related to the sparsity pattern. As the number of wire- wire contacts approaches 0{N'^) in three- 
dimensional packings, the global 6(iV -I- 1) x 6(iV -I- 1) system Jacobian will transform from sparse to 
moderately dense, exposing the iterative NR solver routine to a deteriorating slow-down, assuming that 
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changing contact forces are accounted for. 

Fig. 5 illustrates how the integration schemes compare in terms of relative and absolute serial perfor- 
mance using a small spherical three-dimensional cavity of effective size R/r — 10, with E = 10, v — 0.3, 
(0 = 1, r = 1, /i = 2, Win = 0.005, c = 0.1. The shown curves belong to major disjunct parts of the 
simulation: 

• corotation: computation of the corotational transformation matrices Fg, including Ug via Eq. (20), 
starting from the update of the nodal triads T„ 

• fint assembly: computation of the internal force vector from element contributions (23) 

• Kt assembly: computation of the tangent stiffness matrix from element contributions (24), via (22) 

• foxt assembly: computation of the external force vector from cavity contact forces (28) and self- 
contact forces (42) 

• Joxt assembly: computation of the external Jacobian matrix from cavity contact (36) and self- 
contact (43) 

• adaptive NR line search: adaptive Newton-Raphson line search integration, following Section 3.2 

• adaptive PC: adaptive predictor-corrector integration, following Section 3.3 

For maximal performance, operations that are common to multiple parts, such as the corotation, are 
executed only once, and the results are provided to subsequent parts without contributing to their com- 
putation time. This allows us to directly compare the two algorithms, since, e.g., the subroutine for the 
assembly of the corotation matrices is identical. While it eats up most of the computational costs in 
the explicit APC scheme (Fig. 5a), the implicit method (Fig. 5b) is spending a largely dominant portion 
of its CPU time assembling the tangent stiffness matrix. As expected, the assembly costs for contact 
force Jacobians significantly increases with the packing density, while the force vector assembly is almost 
negligible. The total CPU times are juxtaposed in Fig. 5c. For each scheme, the setup which reaches a 
density of ^ = 0.7 the fastest is shown. Numerical stability limits the constant time step to At = 0.15 
in the explicit case, while the convergence requirements are responsible for the rather small At = 0.07 in 
implicit integration. With adaptivity, larger average time steps can be reached: The fastest settings are 
'Hmin = 10~^, ?7max = lO^'^, resulting in At « 0.28 (red solid), and allowing 5 to 40 line search iterations 
per time step, yielding At « 0.14 on average (blue solid). In summary, explicit PC is more than an order 
of magnitude faster than implicit integration. 

(a) Explicit (b) Implicit (c) Comparison 




Figure 5: Comparison of computational performance. Left: Integrated relative CPU times for the major program compo- 
nents in the explicit (a) and implicit (b) scheme. See the main text for a detailed description of the shown curves. Right 
(c): Total integrated CPU times for the explicit (red) and implicit (blue) scheme, with best adaptive time stepping (solid) 
and best constant time stepping (dashed). The straight solid line is proportional to (jP . 



5. Simulation results 

Previously published experiments and simulations of wires subject to spatial confinement share a 
spherically or cylindrically symmetric container. The goal of the series of simulations presented in the 
following is to investigate the energy levels arising purely from geometry, and to demonstrate that spherical 
coils are not naturally assumed in soft confinements. To this end, two different sets of varying ellipsoidal 
cavities were chosen. The first set consists of isoareal ellipsoids with surface area S = 4000r^, at different 
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aspect ratios {Rx'-Ry'-Rz), imitating an inelastic membrane. The second set of ellipsoids uses the same 
aspect ratios as the first, but is isochoric at volume V = 2.5xl0^r'^ as if elastic cavities were surrounded 
by an incompressible fluid. The set of cavities comprises a sphere (1:1:1), an oblate spheroid (1:2:2), a 
prolate spheroid (3:1:1), and a scalene ellipsoid (3:2:1). In all simulations, CR-RBT is used in combination 
with the explicit adaptive time step integrator. The relevant simulation parameters are E = 10, v = 0.35, 
p = 1, r = 1, h = 3, c = 0.1. The wire is inserted at a sufficiently low speed Win, allowing the wire to 
continuously relax to its local energy minimum. Note that even though the insertion is very slow, a fully 
dynamic treatment of the wire is indispensable, as the large deformation response in the coiled bulk may 
locally occur must faster than Win. This is the case at sudden slip-offs or snap- ins, for instance, where a 
large amount of elastic energy is suddenly released. 

Fig. 6 shows the resulting dense packings for the set of isochoric ellipsoids, and the total bending 
energies of both sets are juxtaposed in Fig. 7. The bending energies are rescaled to the dimensionless 
quantity U\,R/EI using the radii of the respective spherical cavities, R = (3V^/47r)^/^. For comparison 
to the flat scalene limit, the bending energy in two-dimensional ellipses (3:0:1, pseudo-thickness 2r) of 
corresponding surface area or "volume" , respectively, is also plotted. As the contributions from twist and 
longitudinal compression to the total strain energy are minor compared to the bending energy in all ten 
cases, Uh is safe to be used for distinguishing the total internal energy levels. 

Quite naturally, the preferred winding direction imposed by the single short principal axis of the 
oblate confining geometry leads to a highly ordered coil at very low energy (Fig. 6c). While the scalene 
ellipsoid likewise manages to maintain low energy at early stages due to a preferred winding direction 
along the path of minimum curvature, emerging figure-eight shapes (Fig. 6b) then raise the bending 
energy to a moderately higher level in the isochoric setup. The equally limited spatial extents of the 
prolate cavity perpendicular to its long principal axis, on the other hand, do not dictate any directional 
preference to the wire other than along the largest axis, leading to complex bulk dynamics at zero friction 
(Fig. 6a). The bending energy accommodated by the resulting high-curvature loops near the tips clearly 
overcompensates the low energy domain near the center, making the prolate spheroid the least favorable 
of the simulated confinements in terms of strain energy. 

Remarkably, the energy contained in a packed spherical cavity is comparable to that in a prolate 
spheroid. A cut through the packed sphere (Fig. 6e) points at strongly disordered bending at high 
densities inside the shell-like ordered packing as known from Stoop et al. (2011). 




Figure 6: The simulated isochoric ellipsoids with V = 2.5 X 10^ r"^, each at packing density <j) = 0.65. (a): Prolate (3:1;1). 
(b): Scalene (3;2:1). (c): Oblate (1:2:2). (d): Sphere (1:1:1). (e): Clipped sphere. The color represents the rescaled bending 
energy density using L = <j)V/A. The wire radius is halved for clearer representation. 



The bending energy in packed spheres or oblate spheroids with Rx < Ry = Rz can be estimated 
analytically following the ideas of Purohit et al. (2003): 



(65) 
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Figure 7: Rescaled bending energy as a function of the packing density for the different eUipsoidal cavities. The inset shows 
the number of contacts between wire elements. The straight line is proportional to (fy^ . Data from 15 runs with moving 
average filter. Standard errors below 3 are omitted. 



where 



3(/) 



2 r 



1/3 



(66) 



The segment distance dg > 2r, that is in general a function of the packing density, measures the average 
separation of adjacent wire center hues in cross sections perpendicular to the coiling strands. It can be 
extracted from simulation snapshots, or calculated a priori with a numerical solver in some cases (Purohit 
et al., 2003). Motivated by the observation that the density dependence of ds is very weak (Stoop et al., 
2011), we determined it by fitting Eq. (65) to the numerical data in Fig. 7, assuming constant dg, with 
excellent agreement in all cases. 



6. Conclusions and outlook 

In the present paper, a detailed description of a finite element model to simulate the packing of elastic 
wires in two- and three-dimensional cavities has been given. Reddy's locking-free simplified third-order 
beam theory has been employed to model the bending defiections of the wire, offering an improvement 
in precision over recent discrete element simulations. Geometric nonlinearity has been introduced to 
the theory by means of the corotational formulation provided by Crisfield. A quaternion representation 
has been used to allow for arbitrary wire orientations. Interactions of the wire with hard ellipsoidal 
cavities and wire-wire contacts were consistently included into the finite element model. By comparison 
to published numerical results, the correctness of the presented finite element model has been verified. 

Two numerical recipes for integrating the hyperbolic equations of motion in time have then been 
described: The constant-average acceleration method was used as an implicit scheme, and a predictor- 
corrector method served as the explicit counterpart. Adaptive time stepping is crucial in both methods. 
For the former, all internal and external forces were supplemented by their Jacobians for the Newton- 
Raphson iterations, which is a somewhat cumbersome task for self-interacting wire elements. The explicit 
integration method performs very well, clearly beating implicit integration for the following reasons: The 
rough potential energy landscape found in densely packed systems with strong body self-interaction 
calls for relatively small time steps in order to allow fast convergence, and the computational costs 
accompanying the construction and corotation of the Jacobian impair the overall performance of the 
implicit scheme. 

A few isoareal and isochoric configurations have been simulated. Different ellipsoidal shapes yield 
different coiling structures of fully elastic, frictionless wires due to asymmetric spatial confinement, the 
most interesting ones being the scalene and prolate ellipsoids. Published numerical studies concentrated 
on strict spherical or cylindrical confinement and didn't investigate other geometries. The suboptimal 
energy state of spherical coils compared to other ellipsoidal shapes, as it was found in this work, hints 
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at a wire preferring non-spherical packing also when inserted into flexible biomaterial antra as they 
are encountered in real-life applications such as endovascular coiling in the surgical treatment of brain 
aneurysms, for instance. 

Two major physical components determining the morphological characteristics of coiling wires are 
still missing: (i) A plasticity model is needed to account for realistic elasto-plastic behavior. It should 
be noted that the present tangent description (24) is valid only if the internal force may be written as 
fint — K(u) u, which will no longer hold. Explicit integration schemes are unaffected by this restriction, 
however, (ii) The self-contact and cavity contact models also need to be extended by a dry friction 
model to capture the effects of static and dynamic friction. Intrinsic curvature as a third crucial wire 
characteristic, on the other hand, is trivial to incorporate into the presented finite elements by setting 
non-zero local intrinsic angular degrees of freedom in each corotated frame. Together with the suppression 
of axial rotation at the cavity entrance, pre-curvature controls the amount of order in the wire morphology 
in packed three-dimensional cavities. Eventually, we expect to find a strong dependence of the coiling 
morphology on the softness of the cavity, when modeled by a deformablc thin shell (Vetter et al., 2012). 
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